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Abstract 



We consider an adaptive signal control problem on signalized network using Lighthill- 



Whitham- Richards (LWR) model (Lighthill and Whitham 1955 Richards 1956), with 



traffic-derived emission side constraints. We seek to tackle this problem using a mixed 
integer mathematical programming approach. Such a problem class, which we call LWR- 



Emission (LWR-E), has been analyzed before to some extent. As pointed out by Bertsimas 



et al. (2011 ), the mere fact of having integer variables is not the most significant challenge 
to computing solutions of an LWR-based mathematical programming model; rather, it 
is the presence of the nonconvex emission-related constraints that render the program 
computationally expensive. 

As we will show in this paper, it is computationally practical to solve the LWR-E prob- 
lem when it is expressed as a mixed integer program. In particular, we propose a mixed 
integer linear program (MILP) for the LWR-E which can be efficiently solved with com- 
mercial software. The proposed MILP is meant to explicit capture vehicle spillback, avoid 
traffic holding, and reformulate emission-related side constraints as convex and tractable 
forms. The primary methodologies employed by this paper are: (1) a link-based kinematic 
wave model ( Han et al. 2012[ ) ; (2) an empirical and statistical observation between aggre- 
gated emission rate and certain macroscopic traffic quantity; and (3) robust optimization 
techniques. 



1 Introduction 



Traffic signal is an essential element to the management of the transportation network. For 
the past several decades, signal control strategies have evolved from ones developed based on 
historical information, often referred to as the fixed timing plan, to the generation of control 
strategies in which the control system is fully responsive. In the latter case, the cycle lengths 
and splits of the signal are determined based on real-time information. Representatives of 



such signal-control systems are OPAC (Gartner, 1983), RHODES ( Mirchandani and Head 



2000), SCAT (Sims and Dobinson, 1980) and SCOOT (Hunt et al.. 1982). 
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1.1 Mathematical programming approach to traffic signal control 

The performance of a traffic signal control system depends on the optimization procedure em- 
bedded therein. We distinguish between two optimization procedures: 1) heuristic approach, 
such as those developed with feedback control, genetic algorithms and fuzzy logic; and 2) 
exact approach, such as those arising from mathematical control theory and mathematical 
programming. Among these exact approaches, the mixed integer programs (MIPs) are of par- 
ticular interest and has been used extensively in the signal control literature. |Improta and 
Cantarella (1984) formulated the traffic signal control problem for a single road junction as 



a mixed binary integer program. Lo 



(1999a) and Lo (1999b) employed the cell transmission 



model (CTM) ( jDaganzo 1994 1995) and casted a signal control problem as mixed integer 
linear program. In these papers, the author addressed time- varying traffic patterns and dy- 



namic timing plan. In Lin and Wang ( 2004 ) , the same formulation based on CTM was applied 



to capture more realistic features of signalized junctions such as the total number of vehicle 
stops and signal preemption in the presence of emergency vehicles. One subtle issue asso- 
ciated with CTM-based mathematical programs is the phenomena known as traffic holding, 
which stem from the linear relaxation of the nonlinear dynamic. Such an action induces the 
unintended holding of vehicles, i.e., a vehicle is held at a cell even though there is capacity 
available downstream for the vehicle to advance. The traffic holding can be avoided by in- 



troducing additional binary variables, see Lo (1999c). However, this approach ends up with 



a significant amount of binary variables and yields the program computationally demanding. 
An alternative way to treat holding problem is to manipulate the objective function such that 
the optimization mechanism enforces the full utilization of available capacities in the network. 
This approach however, strongly depends on specific structure of the problem and the under- 
lying optimization procedure. Specific discussion on traffic holding can be found in Shen et 



al. (2007) 



This paper is concerned with controlling signalized junctions where the dynamics of vehic- 
ular flows are governed by the network extension of the Lighthill-Whitham-Richards model. 



In particular, we employ the link-based kinematic wave model (LKWM) proposed by Han et 



al. (2012). This model describes network dynamics with variables associated to the entrance 



and exit of each link. It employs a Newell- type variational argument (Newell, 1993; Daganzo 



2005) to capture shock waves and vehicle spillback. Analytical properties of this model per- 



taining to solution existence, uniqueness and well-posedness are provided in Han et al. (2012). 
A discrete-time version of the LKWM, known as the link transmission model, was discussed 
Yperman et al. ( 2005 ) . In contrast to the cell-based math programming approaches where 



m 



the variables of interest correspond to each cell and each time interval, the model proposed in 
this paper is link-based, i.e. the variables are associated with each link and each time interval. 
The resulting MILP thus substantially reduces the number of (binary) variables and hence the 
computational effort. In addition, the link-based approach prevents vehicle holding within a 
link without using binary variables. The network model employed in this paper captures key 
phenomena of vehicular flow at junctions such as the formation, propagation and dissipation 
of physical queues, spillback and vehicle turning. It also considers important features of signal 
control such as dynamic timing plan and time- varying flow patterns. 



1.2 Lighthill-Whitham-Richards model 



Following the classical model introduced by Lighthill and Whitham (1955) and Richards 



(1956), we model the traffic dynamics on a link with the following first order partial dif- 



2 



ferential equation (PDE), which describes the spatial-temporal evolution of density and flow 

(1.1) 



|p(t, x) + lf( P (t, x)) 



where p(t, x) : [0, +oo) x [a, b] — > [0, pj) is average vehicle density, f(p) : [0, p> am ] — y [0, C] 
is average flow. p) am [ s j am density, C is flow capacity. The function /(•) articulates a 
density-flow relation and is commonly referred to as the fundamental diagram. 



Classical mathematical results on the first-order hyperbolic equations of the form (1.1) can 
be found in Bressan (2000). For a detailed discussion of numerical schemes for conservation 
laws, we refer the reader to Godunov (1959) and LeVeque (1992). A well-known discrete 



version of the LWR model, the Cell Transmission Model (CTM), was introduced by Daganzo 



(1994, 1995). PDE-based models have been studied extensively also in the context of vehicular 



networks, with 


a 


list of selected references including 


Bretti et al. ( 


2006) 


|Daganzo|fll995 


); 


Herty and Klar ( 


2003);|Holden and Risebro 


(1995 


: Jin 


fl2003P;|Lebacque and Khoshyaran ||1999[ |2002|). 



; Coclite et al 


(2005); 


(2010 


); 


Jin and Zhang 



Let us introduce function N(-, •) : [0, +oo) x [a, b] — >• R, such that 

— N(t, x) = -p(t,x), Q- t N (^ x ) = f{p(t>x)) 
Function N(t, x) is sometimes referred to as the Moskowitz function or the Newell-curves 



1.2) 



It has been studied extensively, for example by Claudel and Bayen (2010); Daganzo (2005) 



Moskowitz (1965); Newell (1993). A well-known property of N(-, •) is that it satisfies the 



following Hamilton-Jacobi equation 



0_ 
dt 



N(t, x)-f 



d_ 

dx 



N{t, x) 



(t, x) £ [0, +oo) x [a, b] 



;i.3) 



1.3 Traffic operation with environmental constraints 



Unlike many existing literature on signal control which primarily focus on minimizing con- 
gestion level and travel time, this paper address environmental sustainability in traffic op- 
eration and management by consider vehicle emission-related side constraints. However, as 



pointed out by a survey paper Szeto et al. (2012) and the literature therein, environmental 



considerations result in nonlinear and nonconvex feasible sets and objective functions in the 
mathematical programming formulation, adding computational challenges to the problem. As 



a result, heuristic methods (Ferrari 1995) have been developed for these problems. Classical 



methods such as the inner penalty technique Yang and Bell (1997) and augmented Lagrangian 



multiplier technique Yang et al. (2010) have also been used. 

This paper presents a novel approach to circumvent the aforementioned computational 
challenge by re-formulating the emission side constraints as linear constraints, using empir- 
ical observations and robust optimization. Our analysis relies on an observed relationship 
between aggregated link emission rates and the link occupancy, when certain velocity-based 
emission models are employed. Such empirical data are obtained through extensive numerical 
simulation. Detailed description of the simulation and synthetic data is presented in Section 

El 

Despite the strong correlation between the aggregated emission rate and certain macro- 
scopic traffic quantities (e.g. link occupancy), there are non- negligible errors associated with 
such a relationship. These errors are treated as modeling uncertainty and can be naturally 
handled by a technique called robust optimization (RO), which is explained in the next intro- 
ductory section. 
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1.4 Robust optimization 



Uncertainty is usually not negligible in modeling. Errors and perturbations to a deterministic 
model can render an optimal solution in the ideal case suboptimal in implementation. A 
natural approach to capture uncertainty is by assuming that uncertain parameters follow 
certain probability distributions and by employing the notions and methodologies in stochastic 
programming. However, such an approach has arguably two limitations: 1) Exact knowledge 
of distributions is often difficult to acquire, and 2) stochastic programming is recognized as 
highly intractable to solve even with linear objective function and linear constraint functions. 
In view of these challenges, we propose to handle uncertainty in the perspective of robust 
optimization. 

A robust optimization is a distribution-free uncertainty set approach that seeks to minimize 
the worst-case cost and/or to remain feasible in the worst scenario. Compared to stochas- 
tic programming, robust optimization makes no assumption on the underlying distribution 
of uncertain parameters. Moreover, it has been shown to work as a powerful approximation 
to stochastic programming and even probabilistic models with significantly reduced compu- 
tational cost (|Ben-Tal and Nemirovski} fl998| [l999| |2000| |Bertsimas et al.[ |2011a|b[ |Bandi 



and Bertismas 



2012 Rikun, 2011). Although solutions to robust optimization problems can 



be relatively conservative, the conservatism is adjustable with the flexibility of choosing un- 



provided by Bertsimas et al. (2011a). 



certainty sets (Bertsimas and Sim, 2004). A comprehensive review of robust optimization is 



In this paper, we will invoke the robust optimization approach to capture the errors arising 
from statistical learning of the emission model. The goal of this approach is to maximize the 
total throughput of the signalized network while keeping the vehicle emission below a desired 
level even in the worst case. 



1.5 Oganization 

The rest of this paper is organized as follows. In section [2j we present the link-based kinematic 
wave mode (Han et al. 2012[ ). Section [3] presents the basic mixed integer linear programming 
formulation for signalized network, without any emission constraints. Section [4] extends the 
MILP obtained from Section [3] to incorporate emission side constraints, using statistical learn- 
ing and robust optimization in a general setting. In Section[5j we present details of a numerical 
simulation that gives rise to an empirical observation of relationship between certain aggre- 
gated quantities. Such an observation provides support to our RO-based approach. Section [6] 
presents a numerical example, which demonstrates and evaluates the proposed methodology. 



2 Link-based Kinematic Wave Model 

In this section, we present a kinematic wave model on networks, with a triangular fundamental 



diagram for each link. Unlike the cell-based models Daganzo ( 1994, 1995 ), the proposed model 



does not require modeling or computation in the interior of the link. For this reason, we call 
it the link-based kinematic wave model. 



2.1 State variables of the system 

Consider the link represented by an interval [a, b], with b — a = L > 0. In the derivation 
of the LKWM, we select flow q(t, x) and regime r(t, x) as the state variables for the link, 
instead of density. It is obvious that a single value of q corresponds to two traffic states: 
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1) the free flow phase (r = 0); and 2) the congested phase (r = 1). Therefore, the pair 
(q(t, x), r(t, x)) £ [0, C] x {0, 1} determines a unique density value. This simple observation 
gives rise to the following map 

V(-) : [0, C] x {0, 1} -> [0, p> am ], (q, r) m- p (2.4) 



2.2 Riemann problem at a junction with one incoming link 

Extension of the kinematic wave model to a network turns out to be subtle; the issues asso- 
ciated therein are 1) a proper definition of a weak entropy solution at a junction of arbitrary 
topology; 2) uniqueness and well-posedness of the entropy solution. The reader is referred to 



Han et al. (2012); Garavello and Piccoli (2006); Jin (2010); Jin and Zhang (2003) for some 



specific discussion. A junction model can be analyzed by considering a Riemann problem, 
which is an initial value problem with constant datum on each incoming and outgoing link. 
Due to space limitation, instead of a comprehensive discussion on various types of Riemann 
problems, we focus on the Riemann problem for a particular junction, that is, the one with 
one incoming link and n>2 outgoing links. This is because we assume that during one signal 
phase, cars from only one incoming link can enter the junction. 

In order to model vehicle turning, we fix a traffic distribution matrix 



A 



("1,2 "1,3 • • • 



where < cry < 1, i = 2, . .. ,n + 1, Y^=2 a M = 1- The coefficients cry determines how 
the traffic from the incoming link I\ distributes in percentages to the outgoing link 7j. For 
simplicity, we assume A is time- independent. Note that there is no substantial difficulty with 
transforming our modeling framework to deal with time-varying distribution matrices; such 
extension, however, requires additional information on route choices, which is beyond the 
scope of this paper. 

The next theorem characterizes the solution to the Riemann problem at junction with one 
incoming link. 

Theorem 2.1. Consider a junction with one incoming link I\ and n > 2 outgoing links 
I2, ■ ■ ■ , I n +i- F° r every initial data 2/1,0, • • • > Vn+1,0 G [0, C] x {0, 1}, there exists a unique 
n + 1-tuple 

yi, y n+ i g [0, C] x {0, 1} 

where y~i = (q~i, fj) ; such that the solutions to the initial-boundary value problems at the junction 



p(0,x) = 







^(2/1,0) 
V>(yi) 



'mp& x ) + mfM t > x )) 

p(0, x) = ^(yi,o) 
p(t, en) = i/)(yi) 







n+1 



is the admissible weak solution to the junction problem in the sense defined in \Coclite et al 



(2005). In addition, we have the following characterization: the boundary states (qi, fj), i 



1, . . . , n + 1 are given by 



qi 



1, 
0, 
1, 



„max 
?2 



..max 



ln+1 

«1,2 ' "1,3 ' ' Ol,n+l 

*/ n,o = l 

if n.,0 = 0, q\ = qi t o 

if H,o = 0, qi < q lfi 



(2.5) 
(2.6) 
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r; 



'0, 

1, 
0, 



if n,o 
if n,o 
if n,o 



n + l 



1, Qi = 

1, ft < 



lift 
lift 



where 



max 
li 



Qi,0 + rifl(Ci - q ifi ), 
Ci + r ifi (qifl -Ci), 



n + l 



(2.7) 



21 



(2.9) 



Remark 2.2. The quantity qV 1 ^ is the maximum flux an incoming (outgoing) link can send 
(receive) - a quantity identified as demand (supply) by Lebacque and Khoshyaran $1999 , 2002). 



The aforementioned map (y^o) 
Riemann solver, 



,n+l 



i )i=l,...,n+l 



is commonly referred to as the 



sec 



Garavello and Piccoli (2006) for a formal discussion. Theorem 2.1 de- 



scribes the Riemann solver using the new state variables q and r. The verification of theorem 

For junctions with arbitrary topology, the Riemann Solvers are not 
Yet, in the case of one incoming link, we are able to express the 



2.1 is straightforward, 
available in closed-form. 
Riemann solver explicitly. 

In the next subsection, we analyze shock formation and propagation within one single link. 
The location of the shock wave is crucial as it determines the regime variable r at the two 
boundaries of the link. 



2.3 Shock formation and propagation within the link 

We focus on solutions generated by assuming an initially empty network, i.e. y% t o = (0, 0). 
The key to our analysis is the location of a so-called separating shock, which divides each link 
into two zones: free flow zone (r = 0), and congested zone (r = 1). We begin with the fact 
that if the network is initially empty, then there can be at most one separating shock on each 
link 

Lemma 2.3. For every link li and any solution yi(t, x) = (q%(t, x), r»(t, x)) with yj(0, x) = 
(0, 0), the following statement holds: 

1. For every t > 0, there exists at most one x*(t) £ (dj, hi) such that ri(t, x*(t)—) < 

n(t, x*(t)+) 



2. For all x € [aj, bj\, 



Ti(t, x) = 0, if x < x*(t) 
n(t, x) = 1, if x > x*(t) 



Proof. See Bretti et al. (2006). 



□ 

According to Lemma 2.3 the separating shock emerges from the downstream boundary 
bi of the link, and propagates towards the interior of the link. The speed of this separating 



shock is given by the Rankine-Hugoniot condition Evans (2010). 



It is clear that as long as the separating shock remains in the interior of the link /j, the 
upstream and downstream boundary conditions do not interact. Thus the exit of a link remains 
in the congested phase; while the entrance remains in the free flow phase. Consequently, the 



Riemann Solver in Theorem 2.1 is expressed entirely with exogenous parameters Ci and ai . 
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On the other hand, if the separating shock reaches either boundary, it becomes a latent shock. 
Two cases may arise. 

i) The the shock reaches the exit, i.e. x*(t) = bi. In this case, the current link is dominated by 
free flow phase. The boundary condition at x = Oj directly influences the boundary condition 
at x = bi, in a way expressed by 



Qi (t, k 



<?i ( t - —, CJj 



(2.10) 



where Lj is the length of the link, ki is the speed of forward wave propagation. See Figure [T] 
for an illustration. 

ii) The shock reaches the entrance, i.e. x*(t) = a%. In this case, the current link is dominated 
by congested phase. The boundary condition at x = bi directly affects the boundary condition 

at x Qi 

qi (t, ai+) = q t (t-^, bi) (2.11) 
where Wi is the speed of backward wave propagation. 



free flew phase 





congested phase 




case i) 



case ii) 



Figure 1: Example of latent separating shocks. Case i): the separating shock reaches the right 
(downstream) boundary. Case ii) : the separating shock reaches the left (upstream) boundary. 



In either case, the Riemann solver involves boundary flows ( 2.10| ), (2.11), which are en- 
dogenous. The next key step towards the link-based flow model is the detection of these 
two extreme cases. This can be done with a variational method called the Lax-Hopf formula 



(Aubin et al. 2008; Daganzo, 2005 2006; Lax, 1957, 1973 Newell 1993). 



2.4 The variational approach for detecting latent shock 

This section provides sufficient and necessary condition for the occurrence of the latent shock. 



The derivation is omitted for brevity, we refer the reader to Han et al. (2012) for a detailed 



discussion. Define for each link Ii the cumulative entering and exiting vehicle numbers 



N ijUp (t) 



qi(s, (li) ds, -^i,down 



(s, bi) ds 



'0 J 

Recall that the separating shock £*(•) : [0, +oo) — > [a», bi] is a continuous curve in the t — x 
domain. The following theorem provides sufficient and necessary conditions for the occurrence 
of the latent shock. 
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Theorem 2.4. Let N { ( 

Jacobi equation 



: [0, +00) x [ai, bi] be the unique viscosity solution to Hamilton- 
1.3) satisfying zero initial condition, upstream boundary condition iV i)Up (-) 



and downstream boundary condition Ni^ own (•)• Then for all t > 0, 



x*(t) = ai 
x*(t) = bi 



■^i,up \t) ^ ^i,down I t 



Wi 



. I ti t 1 1 T 

+ Pi L 



Ni )Up [t — JT ) — ^i,down (t) 
Ki 



(2.12) 
(2.13) 



Remark 2.5. The significance of criteria (2.12)-(2.13) is that the two extreme cases can 
be detected without any computation within the link. This is because N up (-), Nd OW n{ 



are 



determined completely by the boundary flows. Theorem 2.4 is the key ingredient of the LKWM, 



which allows the network model to be solved at the link level. 

Analytical properties of the LKWM pertaining to solution existence, uniqueness and well- 



posedness are provided in Han et al. (2012). 



3 Traffic signal control problem based on the LKWM 

In this section, the signal control problem is formulated with LKWM. We start with a single 
junction, with two incoming links L\, L2, and two outgoing links ^3, I4 (Figure [2b. Each link 



1 



2 



Figure 2: A signalized junction with two incoming links and two outgoing links. 

is represented by a spatial interval [a», bi], i = 1, 2, 3, 4. The fundamental diagrams for each 
link is given by 

f ( „\ _ \hp [0, p*] 

h\P) — \ , jams ^ ( * 3 am l ' _ ' ' ' 

with Ci = ki p* being the flow capacity. We make the following assumptions 
Al The network is initially empty. 

A2 Drivers arriving at the junction distribute on the outgoing roads according to some known 
coefficients: 

'0:1.3 "1,4 N 



A 

v «2,3 «2,4 

where a« denotes the percentage of traffic coming from link Ii that distributes to outgoing 
link Lj. 

In the problem setting, the flows <&(•)> i = 1, 2 entering links L\ and L2 are known. In 
practice, q~i(-) can be measured at the entrance using fixed sensors such as loop-detectors. 



S 



3.1 Continuous-time formulation 



In this section, we formulate the constraints of the system in continuous time. In Section |3.2| 
we will reformulate the system dynamic as linear constraints in discrete time using binary 
variables. Notice that the discussion in this section can be the building block for extension to 
networks with multiple intersections. 

Let us fix the planning horizon [0, T] for some fixed T > 0. Introducing the piecewise- 
constant control variables itj(-) : [0, T] — > {0, 1}, i = 1, 2, with the agreement that Ui{t) = 
if the light is red for link Ii, and Ui{t) = 1 if the light is green for link Ij. It is convenient to 
use the following set of notations. For i = 1, 2, 3, 4. 



Qi( 

Qi( 

h{ 

n( 
tax ^ 

( 

N U p,i ( 
Ndown,i ( 



~max 
li 



the flow of cars entering link Ij, 
the flow of cars exiting link Ii, 

the binary variable that indicates the regime at x = ai+, 

the binary variable that indicates the regime at x = b{—, 

the maximum flow allowed to enter the link Ii, 

the maximum flow allowed to exit the link Ii, 

the cumulative number of cars that have entered link Ii, 

the cumulative number of cars that have exited link Ii , 

the signal control variable for link Ii, 



Theorem 3.1. The dynamics at the junction (Figure^) with signal control can be described 
by the following system of differential algebraic equations (DAE) with binary variables. 



d 



N up ,i(i) = qi(t), 




d 
dt 



Ndown,i{t) — q-iif), 



if N up j(t) > N down 
otherwise 

if N up , i {t-{ i j _ 
otherwise 

n, [t-—)-Q 

WiJ 

+ fi(t) [Ci-qi 



■ III It r 

+ Pi L 



< N, 



downA 



if Ui(t) 
if Ui(t) 




1 



q k (t) = a x k qi(t) + a 2 ,k fa(t), 



k 

ui(t) + u 2 {t) = 1 for all t G [0, T] 



1, 2, 3, 4 
1,2, 3,4 

1, 2, 3, 4 
1, 2 

:3, 4 

1, 2 
3, 4, 



(3.14) 
(3.15) 

(3.16) 

(3.17) 
(3.18) 

(3.19) 

(3.20) 
(3.21) 



Proof. (3.14) is by definition. For i = 1,2,3,4, if the separating shock on link Ii reaches the 
entrance Oj (exit bi), then the regime variable f j = 1 (fj = 0). Then (3.15)-(3.16) follows from 
Theorem 12.41 

The demand function g™ ax (-) for incoming links and the supply function ^" aa: (-) for outgo- 
ing links are given by (2.9): for i = 1,2, if fi(t) = 1, then g™ a:r (t) = Cf, otherwise if ri{t) = 0, 
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then according to (2.10), 



"max 
li 



(t) = qi(t, h-) = qi [t - -r^, Oj 



Qi t 



U 
k. 



This shows (3.18). One can similarly show (3.17) using (2.9) and (2.11). 



For i = 1, 2, if ui{t) = which means the light is red, then the flow allowed through is 



zero, otherwise, it is given by(2.5). This proves (3.19) 



(3.20) follows from the definition of the splitting parameters 014 i = 1, 2, k = 3, 4. (3.20) 



guarantees that at each time, there are one and only one incoming road that has green light. 

□ 



3.2 Discrete-time formulation 

In this section, we present the discrete-time version of the optimization problem in Theorem 



3.1 Let us introduce a few more notations for the convenience of our presentation. Consider 
a uniform time grid 







t° < i 1 



< t 



N 



t j - t? 



-1 



st, j 



, N 



Throughout the rest of this article, we use superscript 'j' to denote the discrete value evaluated 
at time step V . In addition, we let Li/h = A{St, U/wi = A\8t, A{ G N, A^ G N, i = 
1,2,3,4. 

Approximating the numerical integration with rectangular quadratures, we write equality 



(3.15) and (3.16) in discrete time as 



st E $ - 5t Hi + < mo- r?) 



3=0 
fc-A-f 



A? < k < N, i 



1, 2, 3, 4 



8tY,4 + f^ am Li > -Mr'Z + e 

3=0 



st E i- 5t Y,i ^ M 



fc-A 



/ 



Aj < k < N, i 



(3.22) 

1, 2, 3, 4 (3.23) 



st E i~ 5t 114 ^ M(r*-l) + e 
3=0 3=0 



where f*, f* G {0, 1}. M £ M.+ is a sufficiently large number, e G R+ is a sufficiently small 
number. Constraints ( |3.22 ) and (3.23) determines the regime variables associated with the 
two boundaries of each link. Once the flow phases are determined, the demand and supply 



functions (3.17), (3.18) are re-written in discrete time as 



«zf A *- M (1 -rf) < q™ x ' j < +M{1-H) 



qi A ' Mrj < C aXJ < C Af +Mr> 



1, 2, 3, 4 (3.24) 
1, 2, 3, 4 (3.25) 
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Next, let us re- formulate (3.19). Introducing dummy variables Q 2 i 1 — 3 ' — N, such that 



I ~max,-> 

mm < q • . 



—max, r t -max,') 
03 ?4 

«i,3 «i,4 



Then the discrete-time version of (3.19) can be readily written as 

Jo < qj < Mu\ 
lcj • M(u{ 1) < q\ < Ci 



1, 2 



1,2, j = 1,...,N 



(3.26) 



(3.27) 



In order to write ( 3.26| ) as linear constraints, one could write it as three "less or equal" 
statements, which is simple but bear the potential limitation of traffic holding. Instead, one 
may introduce additional binary variables K , rn and real variables fi\ for i = 1, 2, j = 1, . . . , N, 



such that (3.26) can be accurately formulated as 

,4 Mi\ -!) < ft < qT X Va iA 

<i, ■ m //•/ < ci < <r j '- j 
% -mix- 4) < a <- t; 



-max,-) i 
% M,3 
-max,] i 
QA /Oii,4 
~max,j 



1, 2 



(3.28) 



Finally, we have the obvious relations 

q{{t) = a^ k q{{t) + a2,kq J 2 { t ) 

and 



3, 4, 



1,...,N 



1 



3 



N 



(3.29) 
(3.30) 



The proposed MILP formulation of signal control problem is summarized by ( 3.22 )-( 3.25) 



and (3.27)-(3.30). This formulation captures many desirable features of vehicular flow on 



networks such as physical queues, spill back, vehicle turning, and shock formation and prop- 
agation (although not explicitly). The signal control allows time- varying cycle length and 
splits, as well as the utilization of real-time information of traffic flows. 



3.3 Bound the separating shock 

At the end of this section, we discuss an additional linear constraint that ensures that the 
congested phase on each link is bounded. Such condition is closely related to travel delay: if 
the congested phase remain bounded, there will be a reasonable upper bound for the travel 
time of each driver. Let us consider a single link [a, b]. If we wish to bound the separating 
shock within the interval [c, b] for some a < c < b, this means (with the same notation as 
before) that 



q(s, c) ds < 



t-- 



q(s) ds + p> am (b - c) 



(3.31) 



(3.31) follows by applying (2.12) to the interval [c, b\. It is helpful to notice that since [a, c] 



remains in the free flow phase, q(s, c) must be equal to q (s — ^p). Thus the condition for 
the congested phase to remain in [c, b] becomes 



N up | t 



k 



— -^down I t 



W 



+ P> C 



(3.32) 



Condition (3.32) can be easily written as linear constraint in discrete time. If one wishes to 



include (3.32) in the objective function instead of using it as a constraint, he/she may simply 



minimize the difference between the left and right hand sides of (3.32) 
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4 Emission-related side constraints 



In this section, we propose an emission constraint that employs robust optimization to handle 
parameter errors. It is intuitively clear that the total emission amount on an link of interest 
is highly correlated with the total number of vehicles on the same link. Therefore, it is 
natural to approximate a relationship between the number of vehicles on link I{, denoted by 
Ni(t) = Ni tUp (t) — Ni t d own (t), and the aggregated emission rate e« on link Ij. We assume such 
a relationship is expressible by a polynomial with power C. 

c c 

ei(a, Ni(t)) = ^2ai(Ni(t)) 1 = J^K^W*) ~ N itdown (t)) 1 = a T Ni(t) (4.33) 

1=0 1=0 

where a = (ao, at, . . . , ac) T , Ni(t) = ((A^(t))°, . . . , (Ni(t)) c ). With the these notations, 
the emission side constraint is expressed as 

f €i(a, Ni(t))dt < Ei (4.34) 
Jo 

Such a constraint requires that the emission amount is below some critical value Ei for any 
link I{. Notice that such constraints can be easily transformed to address other types of 
environmental considerations including 

• the total emission amount on the entire network is bounded by some given value 

• the differences among total emissions of all links are bounded or minimized 

The last one ensures that no link (or the neighborhood near that link) suffers much more than 



other links in terms of air quality. Such an issue is identified by Benedek and Rilett ( 1998 ) as 

environmental equity. 

4.1 Side constraints considering errors 



In (4.34), approximation is involved in the relationship between the vehicle number and the 
emission amount. In order to ensure the emission constraint is satisfied even with parameter 
errors, the following robust constraint is of our interest, instead: 

f h{a{t), Ni{t))dt < Ei for all a(t) G (4.35) 
Jo 

where by allowing the coefficient a(t) = (ai(t) : < I < C) to be varied over time, we have 

c 

ei(a{t), Ni(t)) = J2ai(t)(Ni(t)) 1 (4-36) 

fc=0 



and r] a: i is specified as a budget-like uncertainty similar to Atamtiirk and Zhang (2007) given 
as: 



ja: Li < ai(t) < U h for all < I < £, Y^f o a i(t)dt < T ^ =1 ^ j (4.37) 



where L := (Li : < I < C) £ R c+1 and U := (U t : < I < C) £ R c+1 are lower and upper 
bounds. 8 G M+ is a scalar that adjusts the conservatism. Such a constraint corresponds to 
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the notion of robust optimization in the sense that the emission amount on link i is restricted 
to an upper bound with any possible realization of the parameter a(t). 

Increasing the value of 9 leads to a less conservative model. In the most conservative case, 



< 9 < 1, the last constraint in (4.51) is out of effect. 



4.2 Discretization and explicit reformulation 



Constraint (4.35) can be time-discretized into the following form 

N C 



^2^2ai tk (N i)k ) l dt < Ei for all a G f\ a ^ (4.38) 



k=0 1=0 

where a := {a^ k ■ < I < C, < k < N) with the uncertainty set discretized as: 



Va,i 



la: Li< ai )k < U h for all < I < C, < k < N, ^ ^ a tyk 5t < ^ l = l 1 \ 
{ k =o i=i J 



(4.39) 



Such a constraint is in fact a semi-infinite constraint with an infinite index set, which is not 
directly computable. The following theorem gives its computable reformulation. 

Theorem 4.1. Assume that the robust program has a nonempty feasible region. If rj a ^ is 



nonempty, the semi-infinite constraint (4-38) is equivalent to the following set of constraints: 



C N c N N 



1=1 k =o 1=1 k =0 k=0 



s.t. d U)k - d 2 ,z,fc + d 3 = {N i>k y5t for all 1 < I < C, < k < N (4.41) 
di,i, k , d 2 ,i, k , d 3 >0 for all !</<£, < k < N (4.42) 



Proof. Constraint (4.38) can be trivially rewritten as: 

N C N 



max J2J2 a hk( N i,k) l 8t+ max Y^^M^kf^ < E u (4.43) 
which is equivalent to 

max £ Y, ai,k( N i,k? St < Ei - U St. (4.44) 



N C N 



k=0 1=1 k=0 
The evaluation of the constraint function involves solving a parametric problem of the form: 

N C 



max > > ai tk (N Lk ) l 5t (4.45) 

a — ' ' 

k=0 1=1 

s.t. a G fj a>i (4.46) 
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where we treat the Nik as a constant parameter. Program (4.45)-(4.46) has a dual problem 
of the following: 



C N C N 

min di,{,jfct/i - d 2 ,z,fci/ + fy 

Z=l fc=0 z=i k=o 

s.i. rf 1)J)fc - d 2 ,z,fc + ^3 = (Ni ik ) l St for all 1 < Z < C, < k < N 
diik, d 2 ik, d 3 >0 for all 1 < I < C, 0<k<N 



(4.47) 

(4.48) 
(4.49) 



Under the assumption of nonempty feasible region and nonempty f} a ,i-> by noticing the com- 



pactness of fj ai i, the primal program (4.45)-(4.46) and the dual program (4.47)-(4.49) have 



finite solutions and zero duality gap. Moreover, by duality, the objective value of any feasi- 



ble solution to the dual problem ( 4.47 )-( 4.49) provides upper bounds to the primal problem 



(4.45)-(4.46). Therefore, if there exist dij^, c?2,z,fc and d% such that (4.48)-(4.49) are satisfied 



then, if N ik satisfy (4.40), then (4.44) is satisfied. 



On the other hand, if there exists iVj & that satisfy (4.44), then the objective value of the 



optimal solution to the parametric problem ( 4.45 )-( 4.46) is bounded above and thus there 
exists tii / fc, ^2,z,fc and d% such that (4.48)-(4.49) are satisfied. Thus, the equivalence of interest 
is proved. □ 



Remark 4.2. The above theorem is based on the discussions in \Ben-Tal and Nemirovski 
[1999ty and Bertsimas et al. (2011a). With such reformulation the original semi-infinite con- 
straint is now computable with standard nonlinear programming techniques. We refer to the 



reformulated program with constraints (4-40)-(4-4^ as ^ e robust counterpart to the original 
robust problem. 



4.3 A special case when the relationship is linear 

In an desirably simple case where a linear correlation between the vehicle number and the 
emission amount is detected, the robust problem can be considerably simplified. In such a 



case, we can reduce (4.35) to 

r-T 

Mt) (Ni(t)) + a {t)]dt 







where 



[ai(t) (N up>i (t) - N down ,i{t)) + a (t)]dt < E h for all («#(*), of (t)) T G r] a>i (4.50) 



Va,i = { (aftt), a T 2 {t)) T : L\ < m (t) <U u l = 0, 1, £a 1 (t)dt < ^ = {N + * )5tUl 



With discretization, we have the following formulation, 

N 



k=0 



k k 

ai,k ( St ^2 4 - 6t ^2 % ) + a °'' 

3=0 3=0 



St < Ei for all (cjq, tti) € fj a 



(4.51) 



(4.52) 
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where ao := (ao,fc ■ < k < N), a% := (ai jt : < < iV). The uncertainty set is given as 

Va,i = | (ao, 01) : Lq < a ,k < U , L\ < ai jk < Ui, for all < k < N, ^ oi ;fc < — "^g^ 1 | 
I fc=o " ) 

(4.53) 

With the above uncertainty set, we have the following reformulation: 

N N 

- L x d 2 , k + {N + e 1)Ul d 3 + (N + l)U 2 St < Ek (4.54) 

k=0 k=0 

k k 

d 3 + di tk -(h,k = St 2 ^q{ -5t 2 Y,q{ for all k = 0,1,..., N (4.55) 

3=0 3=0 

di, k , d 2 , k , d 3 > for all k = 0, 1, JV (4.56) 
where difc, d 2 k, d 3 are dual variables. In such a linear case, the semi-infinite constraint 



(4.50) is reformulated into a set of linear constraints. Thus, the robust problem can be solved 



efficiently with the state-of-the-art MIP solvers. 

5 An empirical relationship between link occupancy and ag- 
gregated emission rate 

In order to apply the robust optimization techniques from Section |4j we propose to investigate 
the correlation between the link occupancy and the aggregated emission rate (AER). Recall 
that for link i, N up> i(t), Nd OW n,i(t) denotes the cumulative entering and exiting curves, respec- 
tively. We let pi(t, x) : [0, T] x [a^, bj\ — > [0, p^ am ] be the unique entropy solution to the conser- 



vation law ( 1.1 ) which is consistent with the boundary conditions given by Ni up (t), Ni^owni^)- 



In addition, let us assume a speed-based emission model of the form 

e t = T(v) (5.57) 

where et (in gram / hour) is the amount of emission per unit time, v (in mile / hour) is the 
vehicle velocity. The aggregated emission rate on the link i at time instance t is the following. 

AER(t) = / p(t, x)T(v(t, x)) dx t G [0, T] (5.58) 

J at 

where v(t, x) is expressed in terms of p(t, x) via the fundamental diagram. 

We are interested to know if there is any correlation between AER(t), the aggregated 
emission rate on a link level, and N(t), the current link occupancy. If such a relationship 
exists, at least empirically, then we can apply the robust optimization techniques from Section 
[4] to obtain a tractable mathematical program. 

To this end, we employ the speed-based emission model TRANSYT-7F, which was used 
by |Benedek and Rilett| fll998[ ); [Penic and Upchurch| ( |1992| >; |Nagurney et al.| ( |2Q10"| ) 



e T = 26.3009 ■ eXp( °-° m) (5.59) 

v 

where v (in mile / hour) is the vehicle speed, e x (in gram per mile) is the amount of emission 
per unit travel distance. It is useful to notice that the amount of emission per unit time, 
denoted by et (in gram / hour), is determined by 

e t = e x -v = 26.3009 • exp(0.009928w) (5.60) 
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5.1 Numerical simulation 



We conduct a sequence of simulations to obtain sufficient data for the study of correlations 
between the quantities of interest. In order to proceed, we consider a single link with the 
following parameters. 

k = 30 (mile / hour), w = 10 (mile / hour), f? am = 400 (vehicle / mile), L = 0.3 (mile). 

This choice of parameters implies a flow capacity of C = 3000 (vehicle / hour) . We use a time 
grid 

= t° < t 1 ... < t N = T t 3 -t 3 ~ l = St, j = 1,...,N 

As usual, the superscript 'j' indicates the discrete value at the j th grid point. We assume 
D 3 (S J ), j = 0, . . . , N are the time- varying demand (supply) at the entrance (exit) of the 
link. In other words, D J is the number of cars that need to enter the link during the j th time 
interval; whereas S 3 is the maximum number of cars that can advance to the downstream link 
in the j th time interval. In addition, we let each D 3 and S 3 be uniformly distributed random 
variables whose ranges will be specified below. 
We partition the link into cells, 

a = xq < x\ < ... < xm = b xi — xi-i = 5x, I = 1,...,M. 



Then we use the Godunov scheme introduced by Godunov (1959) to compute the spatial 
distribution of density p?, i = 0, . . . , M at each time interval [tj, tj+i]. The aggregated 
emission rate (5.58) is then approximated by 

M 

AER j = 8x^2p{T(vI), j = 0,...,N 

where t^'s are easily calculated from the fundamental diagram, given p^'s. Finally, the link 
occupancy at the j th interval is given by N 3 = 5x YliLo Pi- 

We consider four different scenarios of demand D 3 and supply S 3 of the link. Recall 
C = 3000 (vehicle/hour) is the flow capacity of the link. 

1. {D J } N =0 are uniformly distributed in [0, C/2]; {S 3 } N =Q are uniformly distributed in 
[0, C].~ 

2. {D J } N =Q are uniformly distributed in [0, C/2]; {S 3 } N =Q are uniformly distributed in 
[0, C/2}. 



3. {-D-' j-^q are uniformly distributed in [0, C]; {S 3 }j =0 are uniformly distributed in [0, C]. 

4. {-D- 7 } =0 are uniformly distributed in [0, C]; {S 3 }^ =0 are uniformly distributed in [0, C/2]. 

Each scenario contains 2000 samples. The scatter plots of aggregated emission rate AER 
verses link occupancy are shown in Figure [3] - [6j corresponding to scenarios 1-4, respectively. 
The four scenarios are plotted together in Figure [7| 

The choices of demand and supply imply that scenarios 1- 4 are expected to have increased 
congestion levels. In the first scenario, the traffic is almost always in free flow phase, this 
means the travel speed is a constant, due to the triangular fundamental diagram. Thus the 
total emission rate is proportional to the number of cars (see Figure pi). In both the second 
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Figure 3: Scenario 1. 



Figure 4: Scenario 2. 
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Figure 5: Scenario 3. 



Figure 6: Scenario 4. 



and third scenarios, the demand and the supply are comparable. Therefore we observe similar 
patterns in Figure [4] and [5j where the straight line corresponds to the case where the link 
is dominated by the free flow phase; and the more scattered portions imply the presence of 
both free flow phase and congested phase. In the last scenario, the link is dominated by the 
congested phase, yet the data still imply a linear relationship between link occupancy and 
aggregated emission rate. 

With simple linear regression method, we may approximate the relationship between link 
occupancy and AER with an affine function. We can also obtain the lower and upper bounds 



for the coefficients, as shown in (4.53). 



6 Numerical example 

In this section, we consider the network consisting of two signalized junctions as shown in 
Figure [8j Assume that the inflow at the beginning of links I\, 1% and 1% are given. These 
quantities can be measured, for example by loop detectors. We also fix the upper bound Ei 
for each link Jj so that the total emission cannot exceeds such an upper bound. The adaptive 
signal control problem with emission side constraints is then formulated as a mixed integer 



linear program with constraints (3.22)-(3.25), (3.27)-(3.30) and (4.54)-(4.56). 
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Figure 7: Scatter plot of link occupancy verses AER. 
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Figure 8: Test network with two signalized intersections. 



6.1 Numerical setting 

We assume the same fundamental diagram for all the links, that is, for %= 1, . . . , 7, 



ki = 30 mile/hour, 



Wj 



10 mile/hour, /7? am = 400 vehicle/mile, = 3000 vehicle/hour 



The lengths of all links are set equally to be 0.3 miles. We choose a time grid of 30 intervals and 
a time step of 0.005 hour (18 seconds). The flow entering links I\,l2, 13 are chosen to be time- 
varying functions whose value at each time interval is randomly generated between veh/hour 
and 3000 veh/hour. In addition, to ensure the performance of our signal control strategy, we 
include further constraint that the congested phase must never reach the upstream boundary 
of each link, i.e. no spillback occurs. This is done by invoking a special case of constraint 
d3l32b, with c 



a. 



Regarding environmental consideration, we impose the constraints that the total emission 
emanating from each link, Ii, . . . , J4, does not exceed 150 grams. 
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The MBIP is solved with ILOG Cplex 12.1.0, which runs with Intel Xeon X5675 Six-Core 
3.06 GHz processor provided by the Penn State Research Computing and Cyberinfrastructure. 



6.2 Numerical results 



In order to have a clear visualization of the optimal signal strategy and the separating shocks on 
each link, we use the boundary datum q~i, q\, i = 1, 2, 3, 4 obtained from the optimal solution 



to construct solutions to the Hamilton- Jacobi equation (1.3), using the Lax-Hopf formula 



For the H-J equation, the separating shock no longer represents discontinuity, rather, it is 
displayed as a 'kink' (discontinuity in the first derivative). The Moskowitz functions N{t, x) 
for links ^3, 1± are shown in Figure [To] and 12, respectively. The Moskowitz functions for link 
I2, I a viewed from a different angle is shown in Figure 11 and 12 



One can clearly observe, in each figure, two types of characteristics: forward ones and 
backward ones. The mutual boundaries of the two regimes are identified as the location of 
the separating shock waves, which we managed to deal with implicitly using the variational 
method. It is also clear that the congested regions never reaches the left boundaries, indicating 
that vehicle spillbacks never occur. 






Figure 11: Moskowitz function of link 12- Figure 12: Moskowitz function of link I4. 
We show the emission results in Table 1 below. For comparison purpose, we also compute 
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the optimal signal strategy without any emission-related constraints. It can be seen from Table 
1 that pollutants are mainly concentrated on links I2 and ^3 without emission constraints; and 
the emission amount on I2 exceeds 150 grams. On the other hand, by considering emissioin 
constraints and using RO approach, we managed to control the emission amount on all four 
links to below 150 grams, even under the worst case. 





Link 1 


Link 2 


Link 3 


Link 4 


Without emission constraints 


99.9 


161.4 


144.7 


99.5 


With emission constraints and RO 


109.8 


149.5 


134.8 


111.4 



Table 1 : Link-specific total emissions (in gram) . The values are computed according to the worse case 
over all uncertainty sets of the parameters. 



7 Conclusion 

This paper proposes a novel robust optimization approach to address emission-related side 
constraints, in the formulation of a link-based mathematical programming problem for traf- 
fic signal control. We take advantage of an empirical relationship between the aggregated 
emission rate and the link occupancy, and treat the error as modeling uncertainty which is 
then handled by the RO technique. Such a computational apparatus turns the otherwise 
non-convex emission side constraints into explicit and tractable forms. 

Further research is underway to 1) validate the linear or polynomial relationships among 
certain macroscopic traffic quantities, from both statistical and analytical points of view; 2) 
reveal additional correlations between aggregated emission rate and certain traffic quantities 
that may facilitate the computational efficiency and tractability of the math programs; and 
3) develop heuristic algorithms and computational paradigms for real-time and large-scale 
deployment. 
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